% The issue was that I was dividing the declination data listed as arc
% minutes from north by 10 in the intermagnet.m file. Now this is changed.

% Now correct the data

time_array_min = (datenum(1995,1,1,0,0:9468000-1,30));
a = length(time_array_min);


%load the new variation files
nc_fname=  '/Users/manojnair/projects/obs_mag_data/usgs/BRW_1995_2012.nc';
[x_data, y_data, z_data , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);


nc_fname=  '/Users/manojnair/projects/obs_mag_data/usgs/BRW_1995_2012_Secular_Removed.nc';
[x_data_o, y_data_o, z_data_o , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);


% load the existing secular variation removed data
%nc_fname=  '/Users/manojnair/Downloads/BRW_1995_2012.nc'
nc_fname=  '/Users/manojnair/projects/obs_mag_data/usgs/DED_1995_2012_Secular_Removed.nc';
[x_data_s, y_data_s, z_data_s , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 0);

L = isnan(x_data_s);
data_array = x_data_s(~L);
time_array = time_array_min(~L);
b1 = min(time_array):365:max(time_array)+10;
sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
v=ppval(time_array_min,sp);
x_data_s  = x_data_s - v';        

L = isnan(y_data_s);
data_array = y_data_s(~L);
time_array = time_array_min(~L);
b1 = min(time_array):365:max(time_array)+10;
sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
v=ppval(time_array_min,sp);
y_data_s  = y_data_s - v';        

L = isnan(z_data_s);
data_array = z_data_s(~L);
time_array = time_array_min(~L);
b1 = min(time_array):365:max(time_array)+10;
sp=spline(b1,data_array(1:60:end)'/spline(b1,eye(length(b1)),time_array(1:60:end)'));
v=ppval(time_array_min,sp);
z_data_s  = z_data_s - v';        


%% find the index for the last data available for the exisitind securalr
%variation data
L = isnan(x_data_s);
lastdata = find(L == 0,1,'last');

% remove trend from the new data for index > exisitng data
L = isnan(x_data) | ~isnan(x_data_s);
L(1:lastdata) = 1;
%L = isnan(x_data);

data_array = x_data(~L);
time_array = time_array_min(~L);
%remove the unwanted data from x_data
x_data(L) = NaN;

% now fit
sp = robustfit(time_array,data_array);

v= sp(1) + sp(2) * time_array_min;
%plot

% subplot(311);
% 
% plot(time_array_min, x_data_s,'r');
% hold on;
% plot(time_array_min, x_data - v')
% 
% datetick;

x_data_s(~L) = x_data(~L) - v(~L)';
%x_data_s = x_data - v';

% do this for the other components

L = isnan(y_data_s);
lastdata = find(L == 0,1,'last');

% remove trend from the new data for index > exisitng data
L = isnan(y_data) | ~isnan(y_data_s);
L(1:lastdata) = 1;
%L = isnan(y_data);
data_array = y_data(~L);
time_array = time_array_min(~L);
%remove the unwanted data from x_data
y_data(L) = NaN;

% now fit
sp = robustfit(time_array,data_array);

v= sp(1) + sp(2) * time_array_min;
%plot
% 
% subplot(312);
% plot(time_array_min, y_data_s,'r');
% hold on;
% plot(time_array_min, y_data - v');
% datetick;
y_data_s(~L) = y_data(~L) - v(~L)';
%y_data_s = y_data - v';

L = isnan(z_data_s);
lastdata = find(L == 0,1,'last');

% remove trend from the new data for index > exisitng data
L = isnan(z_data) | ~isnan(z_data_s);
L(1:lastdata) = 1;

%L = isnan(z_data);
data_array = z_data(~L);
time_array = time_array_min(~L);
%remove the unwanted data from x_data
z_data(L) = NaN;

% now fit
sp = robustfit(time_array,data_array);

v= sp(1) + sp(2) * time_array_min;
%plot

% subplot(313);
% plot(time_array_min, z_data_s,'r');
% hold on;
% plot(time_array_min, z_data - v');

z_data_s(~L) = z_data(~L) - v(~L)';
%z_data_s = z_data - v';

%% put back the data

nc_fname=  '/Users/manojnair/projects/obs_mag_data/usgs/DED_1995_2012_Secular_Removed.nc';
[dummy, dummy, dummy , X_ID, Y_ID, Z_ID, obj, ncid] = read_geomag_netcdf(nc_fname, 0, a, 1);


x_data_s(isnan(x_data_s)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
y_data_s(isnan(y_data_s)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
z_data_s(isnan(z_data_s)) = 99999.9;    %no data (this will get converted to 999999 by int32(99999.9*10)
%multiplied by 10 to preserve the decimal points
netcdf.putVar(ncid, X_ID, 0, a,  int32 (x_data_s .* 10) );
netcdf.putVar(ncid, Y_ID, 0, a,  int32 (y_data_s .* 10) );
netcdf.putVar(ncid, Z_ID, 0, a,  int32 (z_data_s .* 10) );

    netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'), 'geospatial_lat', single(lat));

    netcdf.putAtt(ncid,netcdf.getConstant('NC_GLOBAL'), 'geospatial_lon', single(long));

netcdf.close(ncid);



